% Script for HCR mask creation to isolate signal from different cell types

%method is to take masked HCR signal that has puncta signal only (no
%background) and use filters and image operations to create mask.


%% load data and iterate over slices, use folder 0. raw single Z as reference for slice names (it has all slices from batch 2 and 3 only)


clear all;
close all;
fclose all;

%% load data from folder 2. Registered images hafl resolution

path1='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\0. raw single Z\';

% get the folder contents
d2 = dir(path1); %using this directory to extract slice names is useful because it only contains bacthes 2 and 3 which are the ones that were processed in Matlab.
% remove all files (isdir property is 0)
dfolders = d2([d2(:).isdir]==1);
% remove '.' and '..'
dfolders = dfolders(~ismember({dfolders(:).name},{'.','..'}));
path2='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\2.registered after downscaling to half res\';



%% define max filter matrix

% Create a logical image of a circle with specified
% diameter, center, and image size.
% First create the image.
radius=7; % for Max filter
radius2=10; % for Gaussian filter
imageSizeX = 3*radius;
imageSizeY = 3*radius;
[columnsInImage, rowsInImage] = meshgrid(1:imageSizeX, 1:imageSizeY);
% Next create the circle in the image.
centerX = ceil(imageSizeX / 2); % in the middle.
centerY = ceil(imageSizeY / 2);
% radius = 16; %based on profiles

circlePixelsLogical = (rowsInImage - centerY).^2 ...
    + (columnsInImage - centerX).^2 <= radius.^2;
% circlePixels is a 2D "logical" array.
% Now, display it.
% figure
% image(circlePixelsLogical) ;
% colormap([0 0 0; 1 1 1]);
% title('Binary image of a circle');
% axis square;

%% Loop over slices
 f=waitbar(0,'starting')
tic
for i=1:numel(dfolders)
   
     fracprogress=(i-1)/numel(dfolders);
      waitbar(fracprogress, f, {['Masking HCR images'];['Remaining time: ' num2str((1-fracprogress)*(toc/fracprogress)/60) ' minutes']} );
            
   
    currentSlice=dfolders(i).name;
 % load registered images
     Nxph1=imread([path2 dfolders(i).name '\'  'registered_' dfolders(i).name '_Nxph1.tiff']); %Matlab format out for tif is .tiff
     Slc6a5=imread([path2 dfolders(i).name '\'  'registered_' dfolders(i).name '_Slc6a5.tiff']); %Matlab format out for tif is .tiff
     Aldh1a3=imread([path2 dfolders(i).name '\'  'registered_' dfolders(i).name '_Aldh1a3.tiff']); %Matlab format out for tif is .tiff



% -----------------------------------------------------------------

% Nxph1 channel
% apply a circular max filter
Nxph1thresh=50;
Nxph1Filt=ordfilt2(Nxph1,sum(circlePixelsLogical(:)),circlePixelsLogical);  
% figure
% imshowpair(Nxph1Filt*20,Nxph1*100);
% title('max filter')
% Nxph1Filt=imgaussfilt(Nxph1,8);
% 
% figure
% ax1=subaxis(1,2,1)
% 
% imshowpair(Nxph1Filt*500,Nxph1*100);
% ax2=subaxis(1,2,2)
% imshowpair(double(Nxph1Filt>Nxph1thresh),Nxph1*100)
% title('after gaussian')
% linkaxes([ax1 ax2], 'xy')

figure
imagesc(double(Nxph1Filt>Nxph1thresh))
colormap([cmap('gray',100,0,50)]);
caxis([0 3]) %make mask gray
axis square
axis off
xlim([8005 8205])
ylim([5392 5394+200]) 

figure

imshow(Nxph1);
caxis([200 800])
colormap([cmap('green',100,0,50)]);
axis square
axis off
xlim([8005 8205])
ylim([5392 5394+200]) 
set(gcf,'position',[2549 -29 560 420])

Nxph1Mask=Nxph1Filt>Nxph1thresh;

% -----------------------------------------------------------------
% Slc6a5 channel
% apply a circular max filter
Slc6a5thresh=25;
Slc6a5Filt=ordfilt2(Slc6a5,sum(circlePixelsLogical(:)),circlePixelsLogical);  
% figure  %to show max filtered image
% imshowpair(Slc6a5Filt*20,Slc6a5*100);
% title('max filter')
Slc6a5Filt=imgaussfilt(Slc6a5,8); % 8 is sigma or std of gaussian.
% 
% 
% figure  %to show gaussian and thresholding
% ax1=subaxis(1,2,1)
% imshowpair(Slc6a5Filt*500,Slc6a5*100);
% ax2=subaxis(1,2,2)
% imshowpair(double(Slc6a5Filt>Slc6a5thresh),Slc6a5*100)
% title('after gaussian')
% linkaxes([ax1 ax2], 'xy')


figure
imagesc(double(Slc6a5Filt>Slc6a5thresh))
colormap([cmap('gray',100,0,50)]);
caxis([0 3]) %make mask gray
axis square
axis off
xlim([8005 8205])
ylim([5392 5394+200]) 

figure

imshow(Slc6a5);
caxis([200 1500])
colormap([cmap('cyan',100,0,50)]);
axis square
axis off
xlim([8005 8205])
ylim([5392 5394+200]) 
set(gcf,'position',[2549 -29 560 420])

slc6a5Mask=Slc6a5Filt>Slc6a5thresh;

% -----------------------------------------------------------------
%Aldh1a3 channel
% apply a circular max filter
Aldh1a3thresh=25;
Aldh1a3Filt=ordfilt2(Aldh1a3,sum(circlePixelsLogical(:)),circlePixelsLogical);  
% figure  %to show max filtered image
% imshowpair(Slc6a5Filt*20,Slc6a5*100);
% title('max filter')
Aldh1a3Filt=imgaussfilt(Aldh1a3,8); % 8 is sigma or std of gaussian.

% 
figure  %to show gaussian and thresholding
ax1=subaxis(1,2,1)
axis square
% imshowpair(Aldh1a3Filt*500,Aldh1a3*100);
imshow(Nxph1);
caxis([200 1500])
colormap([cmap('green',100,0,50)]);
xlim([8005 8205])
ylim([5392 5394+200]) 

ax2=subaxis(1,2,2)
axis square
imshowpair(double(Nxph1Filt>Nxph1thresh),Nxph1*100)
colormap([cmap('gray',100,0,50)]);
caxis([0 3]) %make mask gray
axis square
axis off
xlim([8005 8205])
ylim([5392 5394+200]) 


title('after gaussian')
linkaxes([ax1 ax2], 'xy')






figure(4)
imagesc((Aldh1a3Filt>Aldh1a3thresh))
colormap([cmap('gray',100,0,50)]);
caxis([0 3]) %make mask gray
axis square
axis off 
box off
xlim([8005 8205])
ylim([5392 5394+200]) 
set(gcf,'position',[2549 -29 560 420])

figure(5)

imshow(double(Aldh1a3));
caxis([200 1500])
colormap([cmap('magenta',100,0,50)]);
axis square
axis off
xlim([8005 8205])
ylim([5392 5394+200]) 
set(gcf,'position',[2549 -29 560 420])

mask3=(Aldh1a3Filt>Aldh1a3thresh);
mask4=mask3([8005:8205],[5392:5392+200]);
imwrite(mask3,[pathSavetemp   'Aldh1a3_Mask.tiff'])




Aldh1a3Mask=Aldh1a3Filt>Aldh1a3thresh;

pathSavetemp=['C:\Users\Tom�s\Desktop\masks for figure examples\'];

%% apply masks


%create Cell type masks based on channel masks
CCMask=and(and(Nxph1Mask,Aldh1a3Mask), ~slc6a5Mask);
LugaroMask=and(and(Nxph1Mask,slc6a5Mask),~Aldh1a3Mask);
GlobularMask=and(and(Nxph1Mask,slc6a5Mask),Aldh1a3Mask);
MLI2Mask=and(and(~slc6a5Mask, ~Aldh1a3Mask),Nxph1Mask);
GolgiMask=and(and(~Aldh1a3Mask,~Nxph1Mask),slc6a5Mask);

%mask channels
CC_Slc6a5=Slc6a5.*uint16(CCMask);
CC_Aldh1a3=Aldh1a3.*uint16(CCMask);
CC_Nxph1=Nxph1.*uint16(CCMask);

Lugaro_Slc6a5=Slc6a5.*uint16(LugaroMask);
Lugaro_Aldh1a3=Aldh1a3.*uint16(LugaroMask);
Lugaro_Nxph1=Nxph1.*uint16(LugaroMask);

Globular_Slc6a5=Slc6a5.*uint16(GlobularMask);
Globular_Aldh1a3=Aldh1a3.*uint16(GlobularMask);
Globular_Nxph1=Nxph1.*uint16(GlobularMask);

Golgi_Slc6a5=Slc6a5.*uint16(GolgiMask);
Golgi_Aldh1a3=Aldh1a3.*uint16(GolgiMask);
Golgi_Nxph1=Nxph1.*uint16(GolgiMask);

MLI2_Slc6a5=Slc6a5.*uint16(MLI2Mask);
MLI2_Aldh1a3=Aldh1a3.*uint16(MLI2Mask);
MLI2_Nxph1=Nxph1.*uint16(MLI2Mask);


% tempfolder
pathSave=['C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\3.Channel colocalization\' currentSlice '\']; %processed channels are here use find string or change file names to all same format.

imwrite(CC_Nxph1,[pathSave   'CC_Nxph1_' currentSlice '.tiff'])
imwrite(CC_Slc6a5,[pathSave  'CC_Slc6a5_' currentSlice '.tiff'])
imwrite(CC_Aldh1a3,[pathSave  'CC_Aldh1a3_' currentSlice '.tiff'])

imwrite(Lugaro_Nxph1,[pathSave  'Lugaro_Nxph1_' currentSlice '.tiff'])
imwrite(Lugaro_Slc6a5,[pathSave  'Lugaro_Slc6a5_' currentSlice '.tiff'])
imwrite(Lugaro_Aldh1a3,[pathSave  'Lugaro_Aldh1a3_' currentSlice '.tiff'])

imwrite(Golgi_Nxph1,[pathSave  'Golgi_Nxph1_' currentSlice '.tiff'])
imwrite(Golgi_Slc6a5,[pathSave  'Golgi_Slc6a5_' currentSlice '.tiff'])
imwrite(Golgi_Aldh1a3,[pathSave  'Golgi_Aldh1a3_' currentSlice '.tiff'])

imwrite(Globular_Nxph1,[pathSave  'Globular_Nxph1_' currentSlice '.tiff'])
imwrite(Globular_Slc6a5,[pathSave  'Globular_Slc6a5_' currentSlice '.tiff'])
imwrite(Globular_Aldh1a3,[pathSave  'Globular_Aldh1a3_' currentSlice '.tiff'])

imwrite(MLI2_Nxph1,[pathSave  'MLI2_Nxph1_' currentSlice '.tiff'])
imwrite(MLI2_Slc6a5,[pathSave  'MLI2_Slc6a5_' currentSlice '.tiff'])
imwrite(MLI2_Aldh1a3,[pathSave  'MLI2_Aldh1a3_' currentSlice '.tiff'])


end


waitbar(1, f,'Done!')
 toc







